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ABSTRACT 


The  design  and  implementation  of  a symbolic  input  and  computation 
package  and  its  application  to  the  development  of  several  new  surface 
interpolation  schemes  are  described.  Capabilities  such  as  the  composi- 
tion of  operators  and  symbolic  di fferentiation  have  been  incorporated 
into  the  system.  This  allows,  in  particular,  the  specification  of 
boolean  sum  projectors.  The  new  schemes  which  have  been  implemented 
include  an  interpolant  to  randomly  spaced  data  and  a "shape  operator" 
which  has  quadratic  precision. 
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INTRODUCTION 

Recent  work  in  Computer-Aided  Geometric  Design  can  be  classified 
into  the  following  three  a.eas  [9]:  graphical  input/output,  man- 
machine  communications  and  mathematical  representation.  In  any  ef- 
fort toward  developing  new  mathematical  schemes  of  surface  represen- 
tation, a computer,  with  its  associated  graphics  hardware  and  soft- 
ware, has  become  an  indispensable  tool. 

However,  despite  the  current  sophistication  of  computer  graphics 
technology,  there  is  still  room  for  improvement  at.  the  computing  end 
of  CAGD  research.  Each  time  a new  surface  interpolation  scheme  is 
conceived,  a significant  amount  of  work  in  the  form  of  mathematical 
manipulations,  analysis  and  subsequent  software  implementation  has 
to  be  done,  before  a computer  display  can  be  produced.  A large  amount 
of  this  work  can  be  viewed  as  error-prone,  time-consuming  overhead 
when  one  considers  that  much  of  the  underlying  mathematical  operations 
and  software  implementations  of  different  polynomial  interpolation 
schemes  are  basically  similar  in  nature.  They  mostly  deal  with  eval- 
uating functions,  taking  partial  derivatives,  composing  linear 
operators,  and  the  like. 

Much  of  the  tedious  symbolic  algebraic  manipulations  heretofore 
done  by  hand  could  be  implemented  as  part  of  a software  system  for 
researching  CAGD  techniques.  Such  an  approach  would  especially  be 
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useful  to  CAGD  in  view  of  the  theorems  developed  by  Barnhill  and 
Gregory  [2,3],  Gordon  [11],  Cohen  and  Ricsenleld  [7],  which  concern 
boolean  sum  interpolants  and  their  interpolation  and  precision  pro- 
perties. These  theorems  provide  a framework  in  which  old  or  new 
interpolants  can  be  "combined"  to  form  other  interpolants  which  will 
possess  the  desirable  properties  and  eliminate  some  of  the  short- 
comings of  the  original  interpolants.  One  of  the  more  powerful  re- 
sults of  boolean  sum  interpolation  theory  is  a composition  theorem 
due  to  Barnhill  and  Gregory  [2]  that  states  the  boolean  sum  of  two 
projectors 


P ® Q = P + q-pg 

has  at  least  the  interpolation  properties  of  P and  the  function  pre- 
cision of  Q.  Here  we  recall  that  an  operator  P is  linear  if 

P(af  + bg)  * aP(f)  + bP(g) 

and  idempotent  if 

P(P(f))  = P(f) 

and  that  a projector  is  a linear,  idempotent  operator  [11].  The 
polynomial  precision  of  a projector  P is  the  set  of  polynomials  which 
P will  reproduce. 

Poeppelmeier  [15]  provides  an  illustration  of  this  theorem  by 
combining  Shepard's  projector,  which  interpolates  to  function  and 
first  derivative  data  at  randomly  positioned  points,  but  has  only 
linear  precision  (i.e. , it  only  reproduces  planes  exactly)  and  the 


- 


Barnhill -Gregory  nine-parameter  interpolant,  which  has  cubic  preci- 
sion, to  obtain  an  interpolant  with  the  interpolation  properties  of 
Shepard's  formula,  but  with  cubic  precision.  Barnhill  and  Gregory 
[2,3],  Cohen  and  Riesenfeld  [7]  have  also  shown  how  boolean  sum 
theory  can  be  used  to  produce  interpolants  free  from  compatibility 
constraints,  the  requirements  for  some  interpolants  that  the  data 
given  over  the  boundary  of  a patch  be  continuous  everywhere  on  the 
boundary  and  that  the  mixed  partials  be  equal  at  the  patch  corners. 

One  of  the  earliest  (and  by  now  classic),  boolean  sum  interpol- 
ants is  the  Coons  Patch  [8].  It  can  be  derived  as  follows:  define 
projectors  P1  and  by 

P7F(x,y)  = (l-y)F(x.O)  + yF(x.l)  (1.1) 

P2F(x,y)  = (l-x)F(0,y)  f xF(l,y) 

(see  Figure  1).  The  boolean  sum  of  these  projectors  is  then  the 
bilinearly  blended  Coons  Patch: 

(P1  ® P2)F(x,y)  = (P7  + P2  - P7P2)F  (1.2) 

+ (l-y)F(x.O)  + yF(x.l) 

+ (l-x)F(0,y)  + xF(l ,y) 

- O-y)[(l-x)F(0,0)  + xF ( 1 ,0)] 

- y[(l-x)F(0,l)  + xF(l ,1)3 

It  can  be  readily  seen  that  this  last  expression,  especially  the  part 
involving  the  tensor  product  P-^F,  can  be  obtained  formally  by  a 
process  that  is  pure  symbol  manipulation. 
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In  order  to  motivate  this  work  further,  we  look  at  the  defini- 
tion of  a few  more  boolean  sum  interpolants.  First,  for  the  bicubic 
Coons  Patch,  define  P^  and  P?  by 

plF(x.y)  = *0(y)f(x,0)  + 4>0(y) F(x,l)  + «t1  (y)Fy(x,0)  + ^ (y)F  (x.l ) 

(1.3) 

P2F(x,y)  = <t0(x)F(0,y)  + U*0(x)F(l  ,y)  + ^ (x)Fx(0,y)  + ^ (x)Fx(l  ,y) 
where 


♦o<‘>  = 

(t-l)2(2t+l) 

. ♦j(t)  = (t-l)2t 

*0(t)  . 

t2(-2t+3) 

. ip,(t)  = t2(t-l) 

are  the  cardinal  basis  functions  for  cubic  He  I'm  i te  interpolation  on 
[0,1].  The  boolean  sum  of  P^  and  P^ 


(pi  « P2)F(x,y)  - fy(y)F{x,0)  + <i'0(y)F(x,l)  + (y)Fy(x,0) 

+ ^(y)Fy(x,l) 


(1.5) 


+ yx)F(0,y)  + i!’0(x)F(l,y)  + ^ (x)Fx(0,y)  + t>}  (x)Fx(l  ,y) 


- ^(y)[^0(x)F(0,0)  + U'0(x)F(l,0)  + *1(x)Fx(0,0)  + ^ (x)Fx(l ,0)] 

- *0(y)C.'0(x)F(0,l)  + ♦0(x)F(l,l)  ♦ ^(x)Fx(0,l)  + *1(x)Fx(l,I)3 


• -----  =p. 
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- *1(y)C*0(x)Fy(OfO)  + -!-0(x)ry(l,0)  + *i(x)Fxy(0,0)  + ^(x)FxyO,0)] 

- <l»1(y)[*0(x)Fy(0,l)  + *0(x)Fy(l,l)  + ^(x)Fxy(0,1)  + *x(x)Fxy(l,l)] 

then  yields  the  bicubic  Coons  Patch. 

Now  we  examine  an  interpolant  over  a triangular  domain  of  defini- 
tion, Nielson's  interpolant  [4].  See  Figure  2.  Define 

PjF  = xF(l-y.y)  + yF(x.l-x)  (1.6) 

P2F  = (Q,  ® Q2)F  = F(0,y)  + F(x,o)  - F(0,0) 

where 

QjF  = F(0,y)  and  Q2F  = F(x,0) 

P^F  interpolates  to  F on  and  P2F  interpolates  to  F on  E|  U E2> 

(R1  ® p2)F  = xF(l-y.y)  + yF(x.l-x)  (1.7) 

+ F(0,y)  + F(x,0)  - F(0,0) 

- x[F(0,y)  + F(l-y.O)  - F(0,0)]  j 

- y[F(0,l-x)  + F(x,0)  - F(0,0)] 

interpolates  to  F all  along  the  boundary  of  the  standard  triangle. 

Barnhill  and  Gregory  have  generalized  this  to  interpolate  cross  j 

boundary  derivatives  [3], 

I 
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All  of  the  above  examples  demonstrate  quite  vividly  that  the 
computation  of  boolean  sums  is  a symbol  manipulation  process.  While 
the  above  cases  can  all  be  done  fairly  easily  by  hand,  examples  a- 
bound  where  the  computation  of  boolean  sums  is  no  longer  a painless 
operation.  Poeppelmeier's  work  [15]  is  a case  in  point.  Recall 
that  he  used  the  Barnhi 1 1 -Gregory  interpolant 


1 1 

UF(x.y)  = l VT^)(1'x)lF0,i(x’0)  + . V’i(T^)(1"X)lF0,i(x,1‘x) 


i=0 


1 y i 

” 

Vp.f" 

+ t *i(Try)0-y) 

i=0  y 

_F(>0(0,y)  - 

c 

. ax1 . 

(O.y) 

_ 

(1.8) 


[(^(f))(0.0)  - (&<§»(0,0)] 


where 


(P2F)(0,y)  = *Q(y)F(0,0)  + 4'1(y)F0J(0,0)  + *Q(y)F(0,l)  (1.9) 

"F  (y ) Fg  -j  ( o » 1 ) 


rsp2Fi 

[“3xJ(0’y)  = y[*J(y)F(o.o)  + *](y)FOJ(o,0)  + *J(y)F(0J) 

(l.io) 


+ ^(y)FOJ(o,i)]  + ^0(y)F1j0(0,0)  4-  ^(y) 


F0.l(°^ 


•py^] 


x=0 


M.0(y)[F,_0(t),,)  - F0i,(0.))] 


M 


[ 


+ (y ) 


r 9Fn  ,(x,l-x)  ■ 

‘HS; — L 


x=0 


The  <^(t)  and  ^'.(t)  are  the  cardinal  basis  functions  for  Hennite 


two  point  Taylor  interpolation  on  (0,1)  given  in  (1.4).  U is  itself 
a boolean  sum  of  2 projectors.  It  was  then  discretized,  i.e., 
changed  into  a form  which  will  accept  discrete  data,  to  obtain 


UF(x,y)  = ( T -o< ) [ '^Q  ( x ) F ( 0 > 0 ) + ^ (x)F1  ^(0,0)  + 7q(x)F(1,0) 


+ <l»1(x)FKO(l,0)] 


+ ^(*O-x)[(l-x)F0j(0,0)  + xFQ  j (1 ,0)  ] 


+ *b(^)[*0(x)F(O,l)  + ^(x)[Fuo(0,I)  - F0J(0,1)] 


+ *0(x)F(l,0) 


+ ^(x)[Flt0(i,o)  - Fq j (1  *0)33 


+ (T^)(1-x)[J<(l-x)CF1  i0(°.l ) + F0J(0,1)] 


+ x[Flt0(l,0)  + Fq j (1,0)] 


4>q(x)F(0,1  ) - *j(x)CFli0(0,l)  - F0J(0,1)]  - j'g ( X ) F ( 1 , 0 ) 


i 1 


*j(x)CFlf0(I.O)  - F0J(1,0)]}] 


+ VT^y)(1’y){(1"y)Fi,o(0’0)  + yFi,o(OJ)  " y^o(y)F(0>0) 


+ 4»j(y)F0>1(o,o) 


+ ^(y)F(o.i)  + ♦|(y)F0tl(o,i)]  - *0(y)F1>0(o,o) 


• *i(y)c'Fo,i(0’0)  ' Fo,i(0,0)  + Fo,i(1’0):i 


' ¥y)EFi,o{CM)  ' F0,l(0’1)] 

*,(y)[-F0i,(0.1)  * ^-F0jl(0,l)  - Fli0(0,l)  * F ,(1,0) 


L 


+ F1>0(1,0)  + 6F(0,1 ) 


+ 4(FU0(0,1)  - F0J(0,1))  - 6F(1,0)  + 2(Flt0(1,0) 


- F0J(1,0))]]} 

- [-F,,0(0,0)  ♦ F1j0(O,l)  ♦ F0i,(0,0)  - F0i,(l,0)]. 


At  this  point  it  should  be  obvious  that  generating  formulas  like  the 
preceding  by  hand  is  a time-consuming,  painstaking  procedure.  In- 
creasing the  complexity  even  further,  Poeppelmeier  then  took  the 
boolean  sum  of  Shepard 1 s projector  S withtheabove  interpolant. 


where 


SF  = A.(x,y)[F(xi,y.)  + (x-x . )F^{ x. ,yi ) + (y-y^ )Fy(x . ,yi )] 


A^x.y) 


III 

j;f«/ 4 lyV) 


e n o 9 

k=0  ft=0  ((x-x,r  + (y-y.r) 

Jt/k 


0.11) 


(1.12) 


to  obtain  a new  interpolant. 

As  new  interpolants  are  constantly  being  developed  which  invar- 
iably require  the  use  of  boolean  sums  at  some  stage,  for  considera- 
tions such  as  compatibility,  interpolation  and  precision,  we  find  an 
increasing  need  to  incorporate  an  algebraic  manipulation  capability 
into  a function  and  interpolant  display  system. 

Toward  the  above  goals,  this  research  has  been  concerned  with 
translating  some  of  the  mathematical  objects  and  operations  required 
in  developing  new  interpolants  into  programming  constructs.  The 
implementation  includes  a command  language  processor,  scanner, 
parser,  symbolic  computation,  and  formula  evaluation  routines.  This 
package  has  been  incorporated  as  a subsystem  to  SURFED,  an  interact- 
ive interpolant  display  and  manipulation  system  implemented  by  the 
University  of  Utah  CAGD  Group.  The  entire  system  runs  on  an  ESS 
Picture  System  connected  to  a PDP  11/45,  taking  27K  words  of  memory . 
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GENERAL  FUNCTIONAL  SPECIFICATION  AND  STRUCTURE 
OF  the  SYSTEM 

Before  we  proceed  to  describe  the  mechanisms  used  in  the  imple- 
mentation of  the  system,  we  should  be  more  precise  about  the  prob- 
lems we  propose  to  solve.  A good  way  to  give  this  description  is 
via  a high-level  functional  specification  of  the  system.  Some  cap- 
abilities in  the  following  list,  particularly  those  that  pertain 
to  interpolant  display  and  manipulation,  and  domain  specification, 
are  part  of  the  original  SURFED  system,  and  should  be  credited  to 
the  work  of  Riesenfeld,  Little,  Herron  and  Dube. 

System  Requirements 

1.  The  system  should  run  in  interactive  mode,  allowing  the  user  to 
enter  interpolants  and  data  from  the  terminal.  The  Picture 
System  display  and  tablet  allow  him  to  specify  different  views 
of  surfaces  displayed.  SURFED  must  also  allow  the  user  to 
Interactively  manipulate  the  shape  of  the  surfaces  being  dis- 
played by  changing  the  data  being  approximated. 

2.  It  should  be  possible  to  enter  an  interpolant  in  symbolic  form 
and  have  it  displayed  as  a surface  in  either  parametric  or 
explicit  form.  The  forms  of  interpolants  accepted  are  rational 
bivariate  polynomials  which  contain  point  functional  data. 
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Partial  and  mixed  partial  derivatives  may  be  necessary  for  des- 
cribing a functional  expression. 

3.  Surfaces  are  defined  for  the  standard  domains  below: 

a.  Rectangular  grid 

b.  triangulated  grid 

c.  Randomly  spaced  data  points. 

4.  Function  and  interpolant  definitions  are  stored  by  the  system 
and  can  be  referred  to  by  name  when  used  in  subsequent  defini- 
tions. This,  along  with  a composition  operator,  allows  the 
definition  of  boolean  sum  projectors.  Thus  expressions  that 
may  appear  repeatedly  in  subsequent  definitions  need  only  be 
defined  once. 

The  block  diagram  of  Figure  3 further  illustrates  the  structure 
of  the  system.  The  symbol  manipulation  capabilities  are  additions 
to  the  "Human  Interface"  and  "Interpolant  Computations"  parts  of  the 
original  SURFED  system. 

A more  detailed  explanation  of  the  way  the  system  functions  can 
only  be  given  by  a user's  manual  of  the  system,  therefore,  we  devote 
the  rest  of  this  chapter  to  such  a manual. 

User's  Manual  for  the  System 

Input  to  the  system  is  in  the  form  of  a command  file.  The 
commands  accepted  are  the  following: 

DEFINE  CONSTANT 
DEFINE  FUNCTION 


DEFINE  PROJECTOR 
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OELETE 

DISPLAY 

EXIT. 

In  addition,  a REDEFINE  can  be  substituted  for  the  above  occurrences 
of  DEFINE. 

We  now  give  a more  detailed  description  of  the  commands.  The 
items  enclosed  in  angle  brackets:  < > are  described  two  sections 
later,  where  we  give  a BNF  description  of  the  expression  syntax  used 
in  the  system. 

Commands 

1.  DEFINE  CONSTANT 

This  command  begins  a section  of  definitions,  each  of  which  must 
begin  on  a new  line  and  end  in  a semicolon.  Each  constant  defin- 
ition has  the  syntax 

<name>  = <integer>  ; 

2.  DEFINE  FUNCTION 

This  begins  a section  of  function  definitions,  each  of  which 
must  be  preceded  by  one  of  the  following  lines: 

GLOBAL 

LOCAL  CIRCULAR 

LOCAL  RECTANGULAR  <integer> 

where  the  integer  is  between  1 and  4.  The  two  LOCAL  definition 
types  indicate  that  the  function  definition  that  follows  is 
identically  zero  outside  either  a circular  or  rectangular  domain. 
This  is  used  in  defining,  for  example,  haystack  functions  for 
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randomly  spaced  data.  For  each  data  point  one  local  circular 
domain,  and  up  to  four  local  rectangular  domains,  may  be  speci 
fied.  Each  function  has  the  form: 

<name>  = <expr>  ; 

DEFINE  PROJECTOR 

This  command  begins  a section  of  projector  definitions,  each 
with  the  form: 

<name>  = <expr>  ; 


REDEFINE  CONSTANT 
REDEFINE  FUNCTION 
REDEFINE  PROJECTOR 

See  Items  1,  2 and  3 for  the  respective  syntax  of  these  commands. 
Note  that  a <name>  may  appear  in  a REDEFINE  definition  section 
only  if  it  has  been  previously  defined. 

DELETE 

This  deletes  all  definitions  from  the  system  except  for  the  last 
one.  This  command  is  used  to  provide  more  free  storage. 

DISPLAY  <name>  <integer> 

This  command  causes  the  function  or  interpolant  named  to  be  dis- 
played on  the  CRT  and  passes  control  from  the  symbol  manipulation 
subsystem  to  the  rest  of  SURFED. 

The  name  given  must  not  have  been  defined  to  be  a constant.  The 
integer  may  be  1,  3 or  4,  depending  on  whether  the  <name>  is  of 
an  interpolant  or  function  that  is  defined  over  randomly  spaced 
data,  triangular  or  square  patches,  respectively. 
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This  command  must  be  followed  by  a line  containing  either  the 
word  PARAMETRIC  or  EXPLICIT,  which  specifies  the  mode  of  display. 

7.  EXIT 

This  causes  a normal  termination  of  the  SURFED  system. 

Command  Fi 1 e Format 

A.  Each  command  must  appear  on  a separate  line.  Blank  lines  are 
allowed.  Only  the  first  three  letters  of  any  command  word  needs 
to  be  specified.  Thus,  DEF  FUN  is  equivalent  to  DEFINE  FUNCTION 
and  LOC  REC  to  LOCAL  RECTANGULAR. 

B.  Each  definition  section  may  contain  an  arbitrary  number  of  defi- 
nitions. 

C.  Definition  sections  may  appear  in  any  order  and  each  type  of 
definition  section  may  appear  more  than  once. 

D.  Constant,  function  and  projector  names  appearing  in  a definition 
must  be  defined  previous  to  the  current  definition. 


Rules  of  Expressions  Syntax 

The  syntax  of  arithmetic  expression  follows  those  of  FORTRAN, 
with  the  following  exceptions: 

1.  Continuation  lines  need  not  be  denoted. 

2.  <expr>:  = <a  FORTRAN  expression  whose  operand  syntax  has 
been  changed  as  in  item  3 below  >|<iter>  <expr> 
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3.  Operand  Syntax 

<operand>:  = <integer>|<name>|<var>|<functional> 
<integer>:  = <digit>|<digit><integer> 

<digi t>:  = 0 j 1 12 | . . . . |9 
<var>:  = U J V 

<name>:  = <letterl>|<letter>  <digit> 

<letter>:  = A|B| . . . |Y|Z 

<letterl>:  = A!...  E|G| . . . |Q|S|T|W| . .. |Z 

<sub-var>:  = <var>(<sub>) | R ( <sub> ) 

<sub>:  = arithmetic  expression  without  parentheses,  and 
whose  operands  may  only  be  <name>  or  <integer>> 
<functional >:  = F(<argl >,<arg2>) |<diff>Ff<argl>,<arq2>) 

• <argl>:  = 0|1 |U(<sub>) 

<arg2>:  = 0 | 1 |V(<sub>) 

<di ff>:  = DU | DV | DUDV | DVDU 
Examples: 


_9 

3u 


c/n  l \ 

I \u,  I / 


niirM  l > 

U J \ U , I / < 


3^7  F(u..Vl)  as:  DUDVF(U(I),  V(I)). 
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4.  Syntax  of  iterated  product  or  sun  notation: 
<iter>:  = '■'op>[<range>]  |<op>[<range'><cond>] 
<op>:  = 

<range>:  = <nan;e>  = <const>,  <const> 

<const>:  = <integer>|<name> 

<cond>:  = <name>^  = <const> 

Thus,  the  following  mathematical  expression: 


n 

E 

i=l 


tvV  / <ui  - vi>  ^7  F(u 

j*l 


j-r 


Vi))l 


would  appear  in  our  command  file  as: 


@[I  = l,N][(U«:i+l)+V(I))/-i[J=l,10,Jv=I] 

[(U(I)-V(I) )+DVF(U( J-l ) ,V( J-l ) ) ]] 

5.  In  addition,  the  composition  operator  CM(<projector  name>, 
<projector  name>)  allows  the  composition  of  two  projectors 
to  be  defined.  Note  that  the  second  projector  in  a composi- 
tion is  currently  restricted  to  not  contain  any  iterated 
product  or  sum  expression. 


Software  Modules  Structure 


This  section  describes  the  major  routines  that  comprise  the 
system.  Figure  4 contains  an  overlay  tree  diagram  of  these  software 


7 


modules. 

The  Command  File  Processor  reads  the  command  file  and  interprets 
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each  command.  When  a definition  is  specified,  it  enters  the  name 
and  type  (e.g.,  PROJECTOR,  LOCAL  CIRCULAR  FUNCTION,  etc.)  of  the 
definition  into  the  symbol  table,  then  calls  SCAN  to  perform  lexical 
analysis  on  the  input  string.  The  parser  parses  the  scanner  output 
using  the  operator-precedence  parsing  technique,  modified  to  take 
care  of  subscript  expressions  and  nested  summation  and  product 
expressions. 

If  the  definition  being  processed  contains  a composition  opera- 
tor, PARSE  calls  COMPOS  to  perform  the  composition.  It  is  here 
that  APPLY  is  invoked  to  apply  a point  functional  to  a previously 
defined  operator  expression.  When  the  functional  being  applied 
involves  the  taking  of  partial  derivatives,  DIFF  is  called  to  per- 
form the  symbolic  differentiation.  SUBUV  is  then  invoked  to  substi- 
tute the  arguments  of  the  functional  in  the  second  operator  of  the 
composition.  Some  basic  simplification  is  performed,  e.g.,  checking 
for  a 0 or  1 multiplicand  or  exponent. 

When  the  Command  File  Processor  encounters  a DISPLAY  command, 
it  passes  control  to  the  SURFED  System,  which  eventually  calls 
EVAL  to  evaluate  the  postfix  definition  of  the  projector  specified, 
and  generate  values  for  the  surface  to  be  displayed.  Since  the 
function  values  are  not  bound  until  evaluation  time,  a new  expression 
need  not  be  created  each  time  the  parameters  for  a surface  are  changed 
during  shape  manipulation.  EVAL  was  written  by  simulating  recursion 
in  FORTRAN  in  order  to  be  able  to  evaluate  nested  summation  and 
product  expressions. 

In  addition  to  the  above,  there  are  many  utility  routines  which 


are  called  by  the  above  modules  and  reside  in  the  lower  nodes  of  the 
overlay  tree  of  Figure  4. 

Before  we  discuss  the  system  in  further  detail,  we  consider  the 
question  of  why  we  did  not  use  an  existing  algebraic  manipulation  sys- 
tem such  as  REDUCE  [13]  as  a preprocessor  to  perform  the  symbolic  com- 
putation tasks  of  our  system.  The  answer  is  that  this  method  does  not 
provide  a sufficiently  high  level  of  user  interaction.  "...  an  alge- 
braic manipulation  system  achieves  its  greatest  effectiveness  if  used 
in  a highly  interactive  man-machine  environment.  The  steps  which  the 
user  takes  in  solving  a problem  very  often  depends  upon  the  results 
of  preceding  calculations  and,  therefore,  a complete  "Program"  in  the 
conventional  sense  is  often  impossible  to  write  a priori."  [13]. 

The  system  we  implemented,  which  integrates  symbolic  computation  and 
surface  display  functions,  does  satisfy  the  interaction  requirement. 


PROGRAMMING  AMD  DATA  STRUCTURING  TECHNIQUES 


USED  IN  THE  SYSTEM 

We  now  discuss  the  major  design  choices  involved  in  the  concep- 
tion of  the  system. 

FORTRAN  was  selected  for  use  in  writing  the  system  for  the  com- 
pelling reasons  of  portability,  accessibility,  efficiency  and  com- 
patibility with  the  existing  host  system  S'JRFED.  It  is  well  known 
to  any  system  designer  that  the  data  structures  to  be  used  deserve 
primary  consideration  in  the  design  process.  Since  the  objects  we 
are  dealing  with  are  of  the  general  rational,  bivariate  polynomial 
form,  how  shall  we  represent  them  inside  the  computer? 

One  method  is  to  use  a square  array,  with  indices  based  on  zero. 
The  coefficients  of  a bivariate  polynomial  can  be  stored  in  the 
entries  of  the  array  with  the  row  number  of  the  entry  representing 
the  power  of  x,  and  the  column  number  representing  the  power  of  y. 
This  method  has  the  advantages  of  efficient  conputation--operations 
such  as  addition,  multiplication,  taking  partial  derivatives  and 
evaluation  are  fairly  straightforward  to  implement  and  fast  to 
execute. 

But  this  representation  of  polynomials  has  been  criticized  for 
wasting  storage  space  in  the  case  of  sparse,  high  degree  polynomials, 
see  [6].  Thus  existing  algebraic  manipulation  systems,  such  as  REDUCE 
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[13]  and  ALT  RAN  [14]  reject  this  method  in  favor  of  a 1 ist  representa- 
tion which,  alongwith  storage  management  routines,  does  keep  memory 
requirements  down  for  operations  wi  th  sparse,  high  degree  polynomials, 
particularly  if  they  are  in  more  than  two  variables. 

Although  a list  representation  generally  leads  to  easy-to-use 
recursively  defined  algorithms  and  elegant,  top-down  program  develop- 
ment, there  is  the  ususal  time  and  space  overhead  associated  with 
using  a linked  list  data  structure,  and  its  requisite  storage  manage- 
ment routines.  Time  and  space  considerations  often  lead  to  the  use 
of  a linear  array  data  structure,  usually  with  the  sacrifice  of 
flexibility  and  elegance. 

Moreover,  there  are  several  characteristics  peculiar  to  the  sys- 
tem that  prompted  the  use  of  a linear  array  to  represent  the  mathe- 
matical formulae.  One  such  consideration  is  that  the  operands  in  our 
expressions  may  be  arbitrary  in  length.  They  range  from  simple  vari- 
ables and  integers  to  subscripted  variables  and  mixed  partial 
derivative  functionals,  e.g.,  DUDVF(U( 1+1 ) ,V( 1+1 ) ) . In  order  for  a 
list  structure  to  accommodate  this,  either  variable-sized  nodes  or 
additional  levels  of  indirection  would  be  required.  Pattern  matching 
and  searching,  required  by  operations  such  as  operator  composition 
and  variable  value  substitution,  can  be  easily  performed  on  a 
1 inear  string  array  . 

Considerations  such  as  the  above  led  to  the  final  decision  to 
store  the  mathematical  expressions  as  a linear  string  array,  in  post- 
fix form,  similar  to  the  internal  form  for  expressions  used  in  many 
compilers  and  interpreters  [12].  Each  token  in  the  string  has  a 
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type  and  value  field.  A postfix  formulation  has  the  well-known  ad- 
vantage of  efficient  evaluation.  It  is  also  conveniently  the  output 
of  the  standard  operator-precedence  parsing  technique  and  it  presents 
the  tree  structure  of  an  expression  in  an  easy  to  access  fashion. 
Thus,  with  this  data  representation,  algorithms  which  either  do  or 
do  not  require  knowledge  of  the  tree  structure  of  an  expression  can 
process  an  entire  expression  with  a linear  scan. 

The  most  difficult  operation  that  is  required  to  be  performed 
on  the  mathematical  expressions  in  this  system  is  that  of  composition 
of  two  operators.  This  often  involves  the  taking  of  derivatives. 

In  the  next  section  a new  symbolic  differentiation  algorithm  is 
described,  which  has  been  developed  to  work  on  an  expression  in 
linear  array  postfix  form. 

A Non-recursive  Algorithm  for  Symbolic  Differentiation  Using  a 
Linear  Array  Data  Structure 

In  this  algorithm  two  stacks  of  pointers,  rather  than  recursion, 
are  used  to  treat  nested  subexpressions. 

A.  Oata  Structures:  A linear  array,  SYMBOLS,  two  stacks  A and  B. 

B.  Assume  the  algebraic  expression  to  be  differentiated  is  in 
syntactically  correct  postfix  form  and  stored  on  SYMBOLS  begin- 
ning at  FIRST  and  ending  at  LAST. 

C.  The  expression  consists  of  tokens  which  are  either  operators  or 

operands.  The  operators  belong  to  the  set  R = {+,  *,  /,  t). 

In  the  description  that  follows,  we  restrict  R to  be  {+,  *,  1} 
for  the  sake  of  brevity.  Extension  of  the  algorithm  to  include 

II . 
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- and  / should  be  obvious.  The  operands  in  the  actual  implemen- 
tation may  be  integers,  variables,  subscripted  variables  or 
point  functionals  such  as  F (u^  .v^ ) , Fu(l,0),  etc. 

We  note  here  that,  since  the  system  currently  only  deals  with 
rational  polynomial  forms,  the  parser  only  accepts  integer  expon- 
ents, since  allowing  exponents  to  be  expressions  would  require 
a log  function  for  purposes  of  differentiation. 

D.  The  differentiated  result  will  be  another  postfix  expression 
stored  beginning  at  SYMB0LS(LAST+1 ) . 

E.  In  describing  the  algorithm  we  make  the  simplifying  assumption 
that  each  operator  or  operand  occupies  one  entry  in  the 
SYMBOLS  array. 

The  basic  stragecy  involves  scanning  the  input  string  from  left 
to  right,  since  it  is  in  postfix  form;  this  corresponds  to  a post- 
order traversal  of  the  expression  tree. 

Step  1.  Begin  by  differentiating  the  first  symbol,  which  has  to  be 
an  operand. 

Step  2.  If  the  end  of  input  has  been  reached,  stop. 

Step  3.  The  types  of  the  next  two  symbols  are  used  to  identify 
three  possibilities. 

The  next  two  symbols  may  be: 

Case  1 : operand-operator 


i 


■ 


i 

t 

j 

i 


We  output  the  derivative  of  this  subtree  and  set  the 
current  input  pointer  to  point  to  the  operator.  Go 
back  to  Step  2. 
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Case  2:  operand-operand 

This  means  the  right  subtree  associated  with  the  left 
subtree  just  differentiated  is  not  a simple  operand  but 
a subexpression.  Wo  set  the  current  input  pointer  to 
the  next  operand,  output  two  nulls  on  the  output  string, 
push  a pointer  to  the  subtree  on  the  input  string  which 
has  just  been  differentiated  onto  stack  A and  push  a 
pointer  to  the  two  nulls  on  the  output  string  onto  stack 
B.  Go  back  to  Step  2. 

Case  3:  operator-whatever 

This  means  that  we  have  just  finished  differentiating 
and  outputting  a right  subtree.  Now  we  can  pop  stacks 
A and  B to  obtain  pointers  to  the  corresponding  left  sub- 
tree on  both  the  input  and  output  strings.  This  is  neces- 
sary because  for  the  * or  / operator,  we  need  to  access 
both  the  differentiated  and  undifferentiated  forms  of 
both  its  left  and  right  subtrees.  Go  back  to  Step  2. 

A more  precise  description  of  the  algorithm  is  given  below 
using  the  programming  language  ALGOL.  In  this  description 
we  assume,  for  simplicity,  that  each  operator  and  operand 
occupies  one  entry  in  the  SYMBOLS  array. 

F.  Util ity  routines: 

1.  PUSHA,  POPA,  PUSHB,  POPS  operate  on  stacks  A and  B. 

2.  Operator  and  operand  are  predicate  functions  which  test  the 
type  of  their  arguments. 

3.  Get  and  Put  gets  from  and  puts  to  a token  at  the  specified 
position  on  SYMBOLS  and  increments  the  given  pointer. 

Diff  differentiates  an  operand. 


4. 


I 
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begin 

I : = F I RST : FREE:=  LAST  + 1 ; 
put  (FREE,  diff  (SYMBOLS(I) ) ) : 
wh  i 1 e ]_<LAST  do 
begin 

II : = I+1 ; 12 : = I +2 ; 

{Check  next  two  syjnbol  types.} 

If  operand  ( SYMBOLS ( II ) ) then 
Lf  operator  (SYMBOLS ( 12 ) ) then 
begin 

{Case  1:  The  next  two  symbols  are  operand-opera- 
tor} 

if  SYM80LS( 12)  = ’+'  then 

begin 

put  (FREE,  diff  (SYMBOLS ( IT ) ) ) ; 
put  (FREE,  '+•) 
end 

else  if  SYMBOLS(IZ)  = then 
begin 

{Output  pointer  to  SYMBOLS  ( 1 1 ) . We  are  not  concerned 
with  distinguishing  between  pointer  symbols  and  con- 
stants here,  and  leave  that  to  lower- level  routines.} 
put  (FREE,  SYMBOLS ( IT);  put  (FREE,’*')' 
put  (FREE, diff  (SYMB0LS( II ) ) ) ; 

{output  pointers  to  SYMBOLS  (I  ) } 
put  (FREE,  I): 


put  (TREE, 
end 


put  (FREE,  '+'} 


' * 


i 1 


else  begin 

{SYMBOLS (12)  must  be  exponentiation  1 1 * } 
recall  tliat  v.'C  only  allow  constant  exponents, 
put  ( FREE , SYMBOLS ( II)); 
put  (FREE, ) ' 

{Output  pointer  to  SYMBOLS ( I ) } 

put  (FREE, I ) ; 

put  ( FREE , SYMBOLS ( 1 1 ) ) ; 

put  (FREE,!);  put  (FREE,'-’); 

put  (FREE,  put  (FREE,'*') 

end 

I;=I2 

end 

else  begin 

{Case  2:  the  next  two  symbols  are  operand-operand 

pusltb(FREE);  pusha(I);} 

put  (FREE,  NULL);  put  (FREE,  NULL); 

I ; = 1 1 : 

put  (DIFF(SYMBOLSO))) 
end 

else  begin 

{Case  3.  The  next  one  symbol  is  an  operator  which 
cannot  be  1 (•' } 
popa  (APTR);  popb  (BPTR); 


I 
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If  SYMBOLS  (II } '+' 

then  put  (FREE,  '+' ) 
else  begin 

{SYMBOLS (II)  must  be  '*'} 

put  (BPTR.I);  put  (BPTR.I) 

put  (FREE.APTR);  put  (FREE, 

put  (FREE, '+' ) 

end 

I:  = I1 

end 

end; 

Note  that  in  the  resulting  postfix  expression  there  are  pointers 
to  the  root  of  subtrees  (subexpressions).  These  subexpressions  are 
recovered  by  using  a copy  operation  and  noting  that  in  a binary  tree, 
the  number  of  non-terminal  nodes  (operators)  is  always  one  less  than 
the  number  of  terminal  nodes  (operands).  Routines  that  do  basic 
simplification  on  the  resulting  expressions  have  been  implemented 
to  reduce  space  requirements. 


[ 
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AN  APPLICATION 


We  now  describe  an  application  of  the  system  which  uses  the 
Barnhill -Gregory  composition  theorem  to  combine  an  interpoiant  and 
an  approximant  to  obtain  a C interpoiant  to  randomly  spaced  data  with 
quadratic  precision.  The  symbolic  processing  capabilities  added  to 
SURFED  permitted  rapid  implementation  of  these  schemes,  since  the 
time  required  to  translate  from  mathematical  ideas  to  computer-accept- 
able form  bas  been  reduced  to  a fraction  of  the  manual  algebraic 
manipulation,  programming  and  debugging  time  previously  required. 
o 

Scheme  1.  A CL  Interpoiant  to  Randomly  Positioned  Data 

Assume  we  are  given  (x^,  y.),  i = l,...,n,  arbitrarily  spaced 

2 

in  R , and  at  each  (x..,  y^ ) we  are  given  function  value  and  tangent 
information  Fix^.y^),  F^x^y. ),  F^x^.y. ).  We  define  below  a C 
interpoiant  to  this  data. 

First,  pre-process  the  data  by  finding  for  each  (x-,y.)  a 
circle  C.  with  its  center  at  (x^.y^)  and  radius  , chosen  such  that 
does  not  contain  any  other  data  point  in  its  interior.  We  can, 

for  example,  choose  R.  to  be  the  distance  from  ( x ^ , y ^ ) to  its 
closest  neighbor. 
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The  distance  function  d^ , that  is  0 at  (x^y^ ) and  1 on  3C. , is 
defined  by 


di (x,y) 


(x  - x.)2  + (y  - yi  )2 


We  can  now  define  a mollifying  function  H ■ at  each  (x^,y^): 


/q(df(x,y))  , 

for  (x,y) 

H-(x,y)  = { 

lO 

el sewhere 

where 

q(t)  = (t  - 

T ) 3 ( 3t  - 1) 

is  the  fourth  degree  polynomial  with  the  following  cardinal  proper- 
ties: 

q(0)  = 1, 

and 

q'(0)  = q(l)  = q ' ( 1 ) = q"(l)  = 0 


To  use  these  H.  as  cardinal  functions  in  the  interpolant,  we 
define  the  truncated  Taylor  operator: 

Li  F(x,y)  = F(x.,y.)  + (x-x.)  Fx(x.,y.)  + (y-y^)  Fy ( x^  .y^ ) 

Then  the  projector  P defined  by 

* l H . (x,y)  L.F(x,y) 
i=l  1 1 


PF(x,y) 


I 

I 

I 
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yields  a (r  surface  that  interpol  .s 


to  function  and  first  derivative. 


| 
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Scheme  2.  A Precision  Operator 

We  now  present  a projector  which  can  be  used  as  a "shape"  operator, 
i.e.,  when  it  is  used  as  the  second  projector  in  a boolean  sum,  it 
helps  to  produce  a projector  which  preserves  the  shape  of  the  given 
data.  It  has  quadratic  precision  if  we  have  9 data  points  given  on 
a rectangular  grid,  and  is  easy  to  compute.  It  is  based  on  the  idea 
that  a multivariate  polynomial  is  its  own  Taylor  series  expansion. 

That  is,  a truncated  multivariate  Taylor  series  can  be  thought  of  as 
an  operator  with  polynomial  precision  up  to  the  term  of  truncation. 

We  specialize  this  to  the  case  of  a bivariate  Taylor  expansion  about 
the  point  (a,b)  and  truncated  after  the  terms  of  second  degree.  Thus 
we  define  the  operator  B: 


B F(x,y)  = F(a,b)  + (x-a)  F (a,b)  + (y-b)  F ( a , b ) 

* y 


2 ? F (a,b) 

♦ (x-a)2  -H, + (y-b)2  


(4.1) 


+ (x-a) (y-b)  Fxy(a,b) 


While  this  operator  yields  quadratic  precision,  it  requires  data  that 
is  almost  never  given,  namely,  the  first  and  second  derivatives. 

We  recall  that  (Conte  and  deBoor,  p.  196)  for  apolynomial  of  degree 

k.  pk(x)» 


Pk(k)(x) 

in 


Pk[x0’xl xk] 


(4.2) 
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whore  the  right  hand  term  is  t k-th  divided  difference  taken  over 
the  points  Xq.Xj , . . . ,x^. 

We  also  note  that  for 

f(x)  = x2,  f ' (a)  * f(a'h)  = f[a-h,  a+h]  (.4,3) 


Since  we  are  only  concerned  with  what  effect  our  operator  A will 

2 

have  on  the  polynomials  1,  x,  and  x , (4.2)  and  (4.3)  imply  that  we 
can  replace  the  derivative  terms  in  (4.1)  with  appropriate  divided 
differences  and  retain  the  quadratic  precision  property. 

Thus  we  obtain  a bivariate  quadratic  precision  operator: 

B F(x,y)  = F(a,b)  + (x-a)  ^^^s 

+ (y-b ) F(a.>b^t) 


+ (x.a)2  F(a+s ,b)  - 2F(a.b)  + F(a_-s,b) 
2s2 


+ (y.b)2  F(a ,b+t)  - 2F(a,b)  + F(a,t 


+ (x-a) (y-b) 


F(a+s,b+t)  - F(a+s,b-t)  - F(a-s,b+t)  +f(a-s,b-t) 


Note  that  B only  requires  function  value  information  at  9 points  on 
a rectangle 


(a,b-t) 

and  that  B interpolates  to  the  5 points  (a,b),  (a,b+t),  (a,b-t), 
(a-s,b),  and  (a+s,b). 


Graphical  Studies 


Although  Scheme  1 is  a smooth  interpolant  to  function  value  and 

derivative  information  at  randomly  positioned  points,  it  does  not 

have  any  precision  property.  This  means  that  while  this  interpolant 

looks  acceptable  for  data  values  close  to  0,  (see  Table  I and  Figures  5and 

6),  it  does  not  have  the  shape-preservation  property  for  larger 

magnitude  data  values.  Thus  for  a second  set  of  data,  (Table  II) 

Scheme  1 (Figure  7)  looks  "bumpy." 

To  remedy  this  problem  we  take  the  boolean  sum  of  Scheme  1 with 

the  shape  operator  of  Scheme  2 (Figure  8)  and  obtain  the  resultant 

interpolant  of  Figure  9.  The  only  "bumpiness"  in  this  last  figure 

is  induced  by  a large  variation  in  the  given  data.  Figure  10 

demonstrates  the  quadratic  precision  property  of  the  interpolant  which 

is  a boolean  sum  of  Schemes  1 and  2,  by  applying  it  to  the  function 
2 2 

x + y . Figures  11  through  18  show  how  Scheme  1 and  its  boolean 

2 2 

sum  with  Scheme  2 behaves  for  the  functions  (x  + y )/(l  + x + y). 


Figure  5 


Positions  of  the  data  Doints  in  the  [0,1]  x 
[0,1]  domain  for  the  explicit  surfaces  of 
Figures  6 through  20.  See  Tables  I and  II. 


Figure  6.  Scheme  1.  Explicit.  Data:  Table  I. 


Figure  9.  Scheme  1 boolean  summed  with  Scheme  2. 
Explicit.  Data:  Table  II. 


Figure  10.  Scheme  1 boolean  summed  with  Scheme  2 
Explicit.  Data:  x2  + y2. 
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Figure  12. 


Scheme  1 boolean  summed  with  Scheme  2.  Explicit. 
Data:  (x^  + y)/(l  + x + y).  The  midpoint  of 
each  stick  indicates  interpolation  at  that 
data  point. 


Figure  15.  Scheme  2.  Explicit.  Data:  1/2  sin  xn/2 
cos  y7T/2 . 


Figure  16.  Scheme  1 boolean  summed  with  Scheme  2.  Explicit. 

Data:  1/2  sin  xn/2  cos  yir/2.  The  midpoint  of 
each  stick  indicates  interpolation  at  that 
data  point. 
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1/2  sin  xn/2  cos  yir/2,  and  1/2  sin  xti  cos  yn,  over  the  area  [0,1]  x 
[0.1].  Note  that  all  of  the  above  figures  are  of  interpolants 
defined  explicitly  over  the  domain  [0,1]  x [0,1]. 

The  data  sets  for  the  figures  mentioned  above  all  contain  nine 
points  that  are  regularly  positioned  in  the  domain  space.  This  is 
a constraint  in  using  Scheme  2 as  a shape  operator.  We  could  remedy 
this  problem  by  doing  parametric  interpolation  to  the  data  in  the 
(x,y,z)-space  and  specify  nine  of  our  points  in  the  domain  (u,v)- 
space  to  be  regularly  spaced.  Figure  19  demonstrates  the  feasibil- 
ity of  this  approach  for  the  data  set  given  in  Table  III. 

Conclusion 

We  have  seen  that,  with  the  addition  of  some  basic  symbolic 
computation  capabilities  to  SURFED,  the  correct  implementation  of  a 
large  class  of  interpolation  and  approximation  schemes  now  becomes  a 
simple  and  routine  matter.  Tedious,  error-prone  and  mechanical  tasks, 
such  as  translating  from  mathematical  formulation  to  program,  compo- 
sition of  operators  and  formal  differentiation  are  no  longer  the 
burden  of  the  researcher. 

For  the  schemes  decribed  here,  the  actual  symbol  manipulation  time 
is  on  the  order  of  a few  seconds.  The  time-consuming  operation  has 
been  generating  values  for  a surface  to  be  displayed.  This  involves 
calling  EVAL  hundreds  of  times  for  a random  data  interpolant.  This, 
perhaps,  suggests  the  need  for  more  sophisticated  simplification 
routines.  Of  more  certain  benefit  would  be  to  speed  up  the 
evaluation  routine,  maybe  recoding  it  in  a lower-level  language. 
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This  is  a one-time  effort  which  would  improve  the  system  performance 
for  the  implementation  of  all  future  schemes.  This  certainly  has  an 
advantage  over  the  previous  method  of  FORTRAN  coding,  manually 
optimizing  and  interfacing  with  SURFED  for  each  new  scheme  being 
developed. 

The  symbolic  computation  part  of  SURFED  currently  occupies  about 
27K  of  16-bit  words  on  the  PDP-11/45.  Better  memory  management  and 
use  of  secondary  storage  would  probably  also  improve  overall  per- 
formance of  the  system. 

Other  useful  extensions  of  the  system  involves  the  specification 
of  expressions  of  more  general  forms.  Foremost  among  these  is  perhaps 
allowing  the  definition  of  transfinite  interpolants , including  the 
identity  operator.  This  may  imply  the  need  for  automatic  discreti- 
zation of  transfinite  schemes,  once  again  done  symbolically.  The 
expression  types  used  may  also  be  extended  to  include  as  part  of 
their  definition  functions  for  which  symbolic  differentiation  can  be 
performed  and  which  are  common  FORTRAN  library  functions,  e.g.,  exp, 
log  and  the  trigonometric  functions.  Projectors  containing  nested 
summation  and  product  expressions  can  perhaps  be  specified  as  the 
second  operator  in  a composition.  This  implies  being  able  to  sym- 
bolically differentiate  such  expressions. 

The  features  implemented  in  this  system  enable  the  research 
mathematician  to  study  graphically  mathematical  surfaces  without 
becoming  involved  in  the  distracting  details  of  software  development, 
interfacing  and  implementation.  This  kind  of  application  of  symbolic 
computation  to  numerical  analysis  should  benefit  future  research  in 
the  area  of  computer-aided  geometric  design,  by  allowing  the  user 


of  this  package  the  facility  to  explore  and  experiment  in  a far-rang- 
ing and  unemcumbered  environment  of  interactive  symbolic  computation. 

This  research,  as  it  brings  together  symbolic  processing  and 
curved  surface  definition,  emphasizes  the  importance  of  the  non-numer- 
ical  tasks  involved  in  computer-aided  geometric  design.  Consequently, 
it  may  provide  further  impetus  for  researchers  in  this  area  to  adopt 
programming  languages  like  PASCAL,  which  are  considerably  better 
suited  for  these  kinds  of  mixed  computational  problems. 
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have  been  implemented  include  an  interpolant  to  randomly  spaced  data  and 
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a ^shape  operator^  which  has  quadratic  precision. 
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